{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 376,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "def sampling(probDistribution,size=10000):\n",
    "    decimal=np.log10(size).astype(int)\n",
    "    p_array = np.array(probDistribution)\n",
    "    p_array_arounded = (np.around(p_array,decimal) * size).astype(int)\n",
    "    if(p_array_arounded.sum()!=size):\n",
    "        print(\"It's not a correct probDistribution:{}\".format(p_array_arounded.sum()))\n",
    "#     else:\n",
    "#         print(\"p_array_arounded:{}\".format(p_array_arounded))\n",
    "\n",
    "    rand_array=np.zeros(size).astype(int)\n",
    "    offset=0;\n",
    "#     print(rand_array)\n",
    "    for i in range(len(p_array_arounded)):\n",
    "#         print(\"i:{}\".format(i))\n",
    "#         print(p_array_arounded[i])\n",
    "        rand_array[offset:offset+p_array_arounded[i]]=i+1\n",
    "        offset+=p_array_arounded[i]\n",
    "        \n",
    "#     print(\"rand_array:{}\".format(rand_array))        \n",
    "    sample_id=np.random.randint(0,size)\n",
    "#     print(\"sample_id:{}\".format(sample_id))\n",
    "    return rand_array[sample_id]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 388,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hmmGenerate(A,B,PI,T):\n",
    "        states=np.zeros(T).astype(int)\n",
    "        observes=np.zeros(T).astype(int)\n",
    "        qi=sampling(PI)\n",
    "#         print(\"qi:{}\".format(qi))\n",
    "        for t in range(T):\n",
    "            states[t]=qi\n",
    "            observes[t]=sampling(B[qi-1])\n",
    "            qi=sampling(A[qi-1])\n",
    "        return (states,observes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 378,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3\n"
     ]
    }
   ],
   "source": [
    "# test sampling\n",
    "total=0;\n",
    "for i in range(5):\n",
    "    s1 = sampling([0.5,0.5],10)\n",
    "    total+=s1-1\n",
    "print(total)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 429,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "qi:1\n",
      "S:[1 2 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 2 1 2 1 1 1\n",
      " 1 1 1 2 1 1 1 1 1 1 1 1 1 2 2 1 1 1 2 1 1 1 2 2 2 1 1 1 1 2 2 2 1 2 1 2 1\n",
      " 1 1 1 1 1 2 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1]\n",
      "O:[2 2 1 1 1 2 1 1 1 1 1 2 1 1 1 2 1 1 1 1 2 1 2 1 2 2 1 1 1 1 1 2 1 1 1 1 1\n",
      " 1 1 1 1 1 1 1 2 1 2 2 1 2 2 2 2 1 1 1 1 2 1 2 1 2 1 1 1 1 1 2 1 1 1 1 2 1\n",
      " 1 1 2 1 1 1 1 1 2 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1]\n"
     ]
    }
   ],
   "source": [
    "#test sampling\n",
    "# s1 = sampling([0.4234,0.1245, 0.4645345],100,2)\n",
    "# print(s1)\n",
    "def f(M, O):\n",
    "#     print(M)\n",
    "#     print(O)\n",
    "    result = np.array([M[x-1] for x in O])\n",
    "#     print(result)\n",
    "    return result\n",
    "\n",
    "N=[1,2]\n",
    "M=['正','反']\n",
    "A=np.array([[0.8,0.2],[0.8,0.2]])\n",
    "B=np.array([[0.8,0.2],[0.5,0.5]])\n",
    "PI=[0.8,0.2]\n",
    "T=100\n",
    "\n",
    "(S,O)=hmmGenerate(A,B,PI,T)\n",
    "print(\"S:{}\".format(S))\n",
    "print(\"O:{}\".format(O))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 430,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20\n",
      "[2 1 1 1 1 1 1 1 1 2 1 1 1 2 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1\n",
      " 2 1 2 2 1 2 2 1 1 1 2 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 2 1 1 2 1 1 1 1 1 1 1\n",
      " 1 1 1 1 1 1]\n",
      "[2 2 2 2 2 1 1 2 2 1 2 1 2 1 2 1 1 2 1 2]\n",
      "12\n",
      "[[66 14]\n",
      " [ 8 12]]\n",
      "[[0.825 0.175]\n",
      " [0.4   0.6  ]]\n"
     ]
    }
   ],
   "source": [
    "# print(len(S[S==2]))\n",
    "# print(len(O[S==2]))\n",
    "print(len(O[S==2]))\n",
    "O1=O[S==1]\n",
    "O2=O[S==2]\n",
    "print(O1)\n",
    "print(O2)\n",
    "print(len(O2[O2==2]))\n",
    "a=np.array([[len(O1[O1==1]), len(O1[O1==2])],[len(O2[O2==1]), len(O2[O2==2])]])\n",
    "print(a)\n",
    "b=np.array([[0.,0.],[0.,0.]])\n",
    "b[0]=a[0]/a[0].sum()\n",
    "b[1]=a[1]/a[1].sum()\n",
    "\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 216,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 216,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.log10(100).astype(int)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
